Thermodynamics of two flavor QCD from imaginary chemical potentials 
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We study QCD thermodynamics in presence of two independent imaginary chemical potentials 
coupled to two degenerate flavors of staggered quarks. Analytic continuation is used to determine 
non-linear susceptibilities, to test the Hadron Resonance Gas (HRG) model below the zero density 
critical temperature, Tc, and to determine the average phase factor of the fermion determinant. 
Deviations from HRG predictions, of the order of a few percent, are clearly visible for temperatures 
T > 0.95 Tc. The determination of non-linear susceptibilities, using different interpolating functions 
for analytic continuation, gives consistent results and in agreement with Taylor expansion compu- 
tations, apart from some systematic effects at or right above Tc. Results for the average phase 
factor are compared with the predictions of Chiral Perturbation Theory; below Tc we are able to 
distinguish the contribution of different hadron states, which is positive (i.e. tends to mitigate the 
sign problem) in the case of baryons. 

PACS numbers: 11.15.Ha, 12.38.Aw, 12.38.Mh 



I. INTRODUCTION 

The study of QCD at finite temperature and baryon 
density has increasing phenomenological interest related 
to the physics of heavy ion experiments and compact 
astrophysical objects. The main open questions regard 
the location and nature of phase transitions in the QCD 
phase diagram, as well as the properties of strongly in- 
teracting matter around the transitions. A reliable an- 
swer to these questions requires at treatment of QCD 
at a non-perturbative level: unfortunately lattice QCD 
simulations, which are the only available tool for a non- 
perturbative study of the theory based on first principles, 
are not possible at finite baryon chemical potential, be- 
cause of the well known sign problem: the QCD fermion 
determinant becomes complex and the probability inter- 
pretation of the QCD Euclidean action, necessary for 
standard importance sampling Monte-Carlo, is lost. 

A number of strategies have been developed to par- 
tially circumvent that problem, like reweighting tech- 
niques [H, 01 13] I the use of an imaginary chemical potential 
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either for analytic continuation [4, S, |6|, 0, H, H, [Mini, 
or for reconstructing the canonical partition function 
Q Ell, Taylor expansion techniquesjlpjp, [H, [H, 
and non-relativistic expansions [2l|, l22ll23j. 

The aim of the present work is that of exploiting 
the method of analytic continuation from an imaginary 
chemical potential to study the properties of hadronic 
matter around the deconfinement transition in QCD with 
two light flavors {Nf — 2). As an improvement with re- 
spect to previous studies based on analytic continuation, 
we introduce to independent chemical potentials, fii and 
fi2, coupled to the two different quark flavors. That is 
equivalent to the introduction of two independent chem- 
ical potentials, fiB and /j,/, coupled respectively to the 
baryon and to the isospin charges B and I3. 

Our strategy will be to determine the dependence of 
the free energy on the two chemical potentials, apart 



from constant terms, by measuring its first derivatives 
with respect to /xi and /i2 (quark number densities) for 
imaginary values of the two variables, and by then fit- 
ting them by suitable functions, to be continued within 
proper analyticity domains. 

One of our aims is the study of generalized suscepti- 
bilities with respect to different conserved charges of the 
model (baryonic, isospin). These quantities are of sig- 
nificant phenomenological interest and have been deter- 
mined till now mostly by the Taylor expansion method. 
We shall compare our results with those obtained by pre- 
vious studies and comment on the efficiency and system- 
atic effects of analytic continuation. In the confined re- 
gion, i.e. below the critical temperature Tc, we shall be 
able to perform a high precision test of the Hadron Res- 
onance Gas (HRG) model, leading to the uncover of vio- 
lations close to Tc. Finally, the knowledge of the depen- 
dence of the free energy on the two independent chemical 
potentials will allow us a study of the average phase fac- 
tor, which gives a direct measurement of the severeness 
of the sign problem. 

Our study is made for QCD with two flavors of unim- 
proved staggered quarks and is based on a standard 
RHMC algorithm. The choice of parameters is taken 
from Ref. [l^. The paper is organized as follows: In Sec- 
tion II we describe the model that we have investigated 
as well as the relevant physical observables; we also dis- 
cuss the symmetries of the model, which are important 
for the choice of the free energy interpolating functions 
to be used for analytic continuation. In Section HI we 
report the technical details of our numerical simulations. 
In Section IV we present results obtained below Tc and 
compare them to the predictions of the HRG model. In 
Section V we report results obtained above Tc. In Sec- 
tion VI and VII we discuss results obtained respectively 
for generalized susceptibilities and for the analytic con- 
tinuation of the average phase factor. Finally, in Section 
VII, we draw our conclusions. 
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II. Nf = 2 QCD WITH TWO INDEPENDENT 
CHEMICAL POTENTIALS AND ANALYTIC 
CONTINUATION. 

QCD with two continuum degenerate flavors is de- 
scribed, in the (rooted) staggered fermion discretization 
of the theory, by the foUowing partition function 

Z{T) = y"l?C/e-^^[^l(detAf[t/])i/2 (i) 

where Sq is the discretized pure gauge action (standard 
Wilson plaquette action in our case) and M is the stag- 
gered fermion matrix describing 4 continuum flavors. Pe- 
riodic (antiperiodic) boundary conditions are assumed 
for gauge (fermion) fields along the Euclidean time di- 
rection. 

The introduction of two independent chemical poten- 
tials, /zi and 112, coupled to the number operators of each 
quark family leads to the following expression for the 
grand canonical partition function: 

Z(T,^i,^2) = yX'[/e-^^detM3[//i]detM3[^2] (2) 

where the fermion matrix in the standard staggered for- 
mulation at finite chemical potential reads: 

1 ^ 

+ r^^A (e'^^t/..4<5,,,_4 - ^-"'Ul-^.A^^-^ (3) 

Here i and j refer to lattice sites, i) is a unit vector on 
the lattice, rji^^, are staggered phases; a/i and am are re- 
spectively the chemical potential and the quark mass in 
lattice units. 

The two chemical potentials can be rewritten in terms 
of a quark number chemical potential /i^ = (/ii -f /i2)/2 
(or equivalently a baryon chemical potential /is — 3fiq) 
and of an isospin chemical potential = (/ii — /Z2)/2. 

While the original theory is invariant under both 
charge conjugation and isospin rotations, the theory in 
presence of finite chemical potentials obviously is not. 
However the original invariance is reflected in the fact 
that the free energy F = —T\nZ must be an even func- 
tion of jiq and Hi separately, or equivalently it must 
be invariant under the two following transformations 
(/ii,^2) (Ai2,Mi) and {^.1,^.2) (-Ai2,-A*i), which 
are easily verified to be symmetries of the partition func- 
tion in Eq. That places strong constraints on its 
possible functional dependence. 

In presence of a finite chemical potential det M be- 
comes complex and det Af[—/x] = (detM[/i])*. Therefore, 
apart from the case /^2 — (fJ-q — 0), the integrand in 
Eq. ([2]) is complex and cannot be interpreted as a prob- 
ability distribution over gauge fields, so that standard 
importance sampling techniques cannot be applied (sign 
problem). 



Positivity is recovered if the chemical potentials fii and 
fi2 are taken as purely imaginary: in this case numerical 
simulations are feasible and results can be used to fit the 
functional dependence of relevant observables. 

Due to the above mentioned symmetries of the free 
energy, analytic continuation is actually a continuation 
from negative to positive values of fj,^ and fij. Of course 
it is expected to be applicable as long as no phase tran- 
sitions are met along the continuation path. 

It is convenient for the following discussion to intro- 
duce the variables 

9q = Im(^q)/T = 7VtaIm(/Zq) 
01 = lmUii)/T = Ntalmifii) 

and 

ei = im{^ii)/T ^Oq + ei 

O2 = Ini(/i2)/T = 9q-ei 

where Nt is the number of lattice sites in the temporal 
direction. 

It can be easily shown that the introduction of an imag- 
inary chemical potential is equivalent to a twist in the 
temporal boundary conditions for fermions by an angle 
Im(/i)/r. Hence both determinants appearing in Eq. ([2]) 
are periodic functions, respectively of 9i and 62, with pe- 
riod 2tt, so that the free energy itself is a periodic function 
of these variables. 

In terms of 6q and 9i that means again periodicity 
with period 2tt in both variables, plus invariance under 
{6q, 9]) {6q + Tr,6j + tt). However, following the argu- 
ment given by Roberge and Weiss in Ref. [ISj, it is pos- 
sible to prove that a transformation 9q ^ 9q + 2TTk/Nc, 
where Nc is the number of colors and k is an integer, can 
be cancelled by a change of variables in the functional in- 
tegration in which all temporal links at a given time slice 
get multiplied by a center element exp(— ?2fc7r/7Vc) (cen- 
ter transformation). Hence the free energy is expected 
to be a periodic function of 9q with period 2tt/Nc in- 
stead of 27T {Nc = 3 in our case). An analogous change 
of variables does not work for translations in 9i, which 
rotate the link variables appearing in each determinant 
in a different way, therefore the period in 9j is really 2tt. 

For temperatures below the zero density critical tem- 
perature, Tc, no phase transitions are expected, as in 
the hj — case, in the whole 9q,9j plane. Therefore, 
due to the discussed periodicity and required symmetries, 
the most natural parametrization of the free energy is in 
terms of a trigonometric series as follows: 

Ei^^ ^ J2 WH,i cosiihOq) cos{l9i) (4) 

h,l 

with h and / both integers; moreover h and I must 
have the same parity because of the invariance under 
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{dq,6i) {9q + 7T,6i + ir). Further constraints on the 
number of terms appearing in Eq. ^ may be predicted 
by particular effective models of strong interactions be- 
low Tc, like for instance the HRG model to be discussed 
in Section IIVI In such regime, information valid for an- 
alytic continuation can be gathered in the whole 9q,9i 
plane. 

For T > Tc we expect instead phase transitions in the 
9q,9i plane, corresponding either to the continuation of 
the physical deconfinement transition or to the general- 
ization of Roberge- Weiss (RW) transitions. Therefore a 
limited region around 9q = 9i = is available for the pur- 
pose of analytic continuation to real chemical potentials, 
and we shall write an expression for the free energy valid 
in that region which respects the predicted symmetries 



under 9n 



-9„ and 9i 



-9 J separately. In particular 



the free energy will be expressed as a polynomial like 



Fi9q,9l) 

T 



9f 9f 



(5) 



with i, j non negative integers, or as a ratio of polynomi- 
als of the same kind 



F{9q,9l) 

T 



"00=0 



g2k 



k,l (2k)\ (2l)\ 



(6) 



doo = l 



The latter is an example of Chisholm approximant, i.e. 
the generalization to the case of two independent vari- 
ables of usual Fade approximants, which have revealed 
to be better suited for analytic continuation in some 

cases [i,!!!,!!!. 

Some of the quantities we are interested in are general- 
ized susceptibilities with respect to the different chemical 
potentials, which for Nf = 2 are defined as follows 



■P 



(7) 



where P is the pressure. Analogous susceptibilities are 
defined in terms of /j,q and /i/ 



^^i'q^^l] 



(8) 



The free energy symmetries discussed above imply pre- 
cise constraints on the susceptibilities computed at zero 
chemical potentials. In particular we have xt'j only 
if i and j are both even, while Xi,j ^ if z + j is even 
and = Xj.i- 

Such quantities encode all relevant information about 
fiuctuations of conserved charges, which are generally 
considered to be sensitive probes for the properties of 
the thermal medium produced in heavy ion collisions. 
They have been computed mostly in the Taylor expan- 
sion approach [H, [13, IH, IH, where they are ex- 
pressed as average values at ^ = of operators which are 



more and more complex and computationally demand- 
ing as the order grows, since they require more and more 
matrix inversions. It is therefore sensible to explore the 
consistency and the efficiency of different strategies. In 
the analytic continuation approach we determine numer- 
ically the functional dependence, for imaginary values of 
the chemical potentials, of the first derivatives of the free 
energy. In terms of adimensional quantities, which are 
most conveniently determined on the lattice, they are 
given by 



ni 



(Ni) _ d 



(P/T^ 



{P/T^ 



ni + n2 



ni - n2 



(9) 



where Nq and Nj are respectively the quark number and 
isospin charge operators, with 



1 dlnZ 



1 dF 



Tr(^Af-i[C/,Ai,]^M[C/,M.])^ (10) 
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for i ~ 1,2. In terms of the susceptibilities defined in 
Eq. O fii = Xifi/T^ and 712 = Xos/T^- Such first 
derivatives, which are purely imaginary for imaginary 
chemical potentials, can be measured quite efficiently 
(only one matrix inversion is needed for the noisy es- 
timation of the trace) and, apart from constant terms, 
encode all information about the dependence of the free 
energy on /x^, /x/. Information gathered at imaginary val- 
ues of fiq/i can then be analytically continued to real 
values of Hq/i, in particular higher order derivatives at 
pi-q = l^i = ^ can be extracted. 

In comparison to the Taylor expansion approach, the 
great advantage related to the much simpler observables 
can be compensated by the need for multiple simulations 
at different values of the chemical potentials. Moreover, 
this procedure involves some systematic dependence on 
the function chosen to interpolate data at imaginary /i's, 
which should be eventually checked by comparing re- 
sults obtained with different functions. We shall com- 
pare trigometric expansions with polynomials below Tc, 
polynomials with ratio of polynomials above Tc- 



III. PARAMETER DETAILS AND NUMERICAL 
SETUP 

Since we want to compare our results for the gener- 
alized susceptibilities with those obtained by the Taylor 
expansion approach, we have chosen for this study a sub- 
set of the parameters used in Ref. [l^ , which is reported 
in Table HI That corresponds to five different tempera- 
tures with a standard staggered lattice discretization on 
iVi = 4 lattices and a fixed value (on the correspond- 
ing T = lattices) for the pion mass, m^r — 280 MeV 
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given by 



TABLE I: List of simulated temperatures and corresponding 
Pjiriq values (taken from Ref. [l9|), number of (p,q,^i) pairs 
explored at each temperaure {ripairs) and average number of 
trajectories (of 1 MD time length each) generated at each 
temperature and for each {fiq, ^i) pair {ntraj)- Nn instead 
indicates the total number of Dirac matrix multiplications 
performed at each T, which is reported as an estimate of 
the overall computing effort performed: that is more or less 
equally distributed between Monte-Carlo and measurements. 



(actually m^/nip = 0.31(1) and rUp/Tc = 0.54(2)). The 
critical temperature reported in Ref. [l^ is Tc — 170 
MeV. 

In particular, we have made simulations on a 16^ x 4 
lattice using a RHMC algorithm. Our spatial size Ls = 
16 corresponds to about 6.6 inverse pion masses, hence 
finite size effects are not expected to be important. 

For T < Tc we have made simulations on a grid of 
about 100 different pairs {9q,9i), in the range [0,7r] x 
[0, tt]: because of the above described periodicity, this 
surely contains all possible information available at imag- 
inary chemical potentials (actually in a redundant way, 
which however is a benefit for checking the reliability of 
our statistical analysis). Since susceptibilities are calcu- 
lated at null values of /ig and /i/, more points were taken 
in a restricted region around the origin, in order to per- 
form fits of low-degrees polynomials in Uq and nj around 
the origin easily. Morover, we have decided to perform a 
more accurate study of HRG model along the axis 6*/ = 0, 
therefore we have chosen further points there. 

For T > Tc we have performed a preliminary study 
aimed at finding the position of transition lines, with the 
purpose of delimiting the region at imaginary chemical 
potentials available for analytic continuation. Further 
information about this region are given in Section fVl 

For each (T, 9q,9i) we have produced about 2-3K ther- 
malized trajectories of 1 Molecular Dynamics time length 
each. More details about the amount of (/x,, m) pairs ex- 
plored and average numbers of generated configurations 
are given in Table HI 

Quark densities have been measured by using noisy es- 
timators. It is possible to minimize the total error of 
these observables (sum of statistical and noise fluctua- 
tions) at fixed simulation time by choosing an appropri- 
ate number of random vectors used for each noisy esti- 
mation. Assuming that noise and statistical fluctuations 
are independent of each other, the optimal number of 
random vectors riycc to be used for each configurations is 



noise / '^conf 



(11) 



where Umeas is the variance of the observable (quark den- 
sity) over different configurations, anoise is the variance 
of the different estimates of the observable over a fixed 
configuration, Tconf is the time needed to generate a new 
configuration and TesUm is the time needed to perform 
one noisy estimate of the observable. We have measured 
those quantities in preliminary runs and we have found 
that, with our numerical setup, this number is around 30 
for all explored parameter sets. Notice that Eq. pT|) does 
not take into account the autocorrelation among config- 
urations and thus overestimates n^ec! we have however 
directly checked, by comparing different choices of n^ec, 
that the efficiency is almost stable for riycc ~ 10 — 50. We 
have always chosen n^ec = 30 in our production runs. 

Simulations have been done on two PC farms in Genoa 
and in Bari. The complete collection of our data is not 
reported here, but is at disposal for interested readers. 



IV. RESULTS AT T < T^: PRECISION TEST ON 
THE HADRON RESONANCE GAS MODEL 

The thermal medium below the critical temperature is 
generally believed to be well described as a gas of free 
hadron resonances (HRG model). This model provides a 
good description of thermal conditions at freeze-out [1^, 
Wh . and has received theoretical support from lattice 
QCD simulations [l^. Deviations from the model have 
been recently detected close to Tc in a lattice study based 
on the Taylor expansion method . 

In the HRG model the free energy is expressed as the 
sum of free particle energies. In particular, the free en- 
ergy for species i of spin Si, mass rrii, baryon number Bi 
and isospin I^i, is given by 



T\nZ., = ± 



27r2 



27r2 




g^VT^mj ^ [(±l)'+i 



1=1 



where gi = 2si + 1, the upper (lower) sign applies to 
mesons (baryons) and 



exp 



(13) 



The expression in Eq. (|12p is an approximation in the 
case of unstable particles, for which an integration over 
a Breit-Wigner distribution in the particle mass would 
be more appropriate. The Bessel function K2 is expo- 
nentially suppressed for large values of the argument, 
K2{x) ~ •\/7r/(2a;)e~^, hence for mi ^ T we can keep 
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just the first term I = 1 in the I expansion, correspond- 
ing to the Boltzmann approximation in which quantum 
statistics effects are neglected. Summing up over all 
known particles and resonances and grouping together 
all charge conjugation and isospin partners we get 

liiZ = VT^ W{m,g,T)S{B) cosh (sB^'j 

B,I,m 

^ ^(/3) cosh (2/3^) ] (14) 

/3>0 J 

where S(n) = 1 — l/2(5„.o 

Such prediction is easily continued to imaginary chemical 
potentials, where hyperbolic functions get transformed 
into trigonometric functions, in particular we have 

liiZ = VT'^ J2WBjiT)d{B)cos{W9q) 

B,I 

^ ^(J3)cos(2/30/) (15) 



Im(ng) = "Y3BWBjiT)sm{W9g) 

B,I 

^ ^(/3)cos(2/30/) (16) 



Latter prediction is easily affected by lattice artifacts 
and by the unphysical quark masses used in simulations, 
which change the actual hadron spectrum on the lattice. 
The former, instead, is expected to be more robust and 
less sensitive to discretization details. The method of 
analytic continuation is particularly well suited for lat- 
tice QCD tests of the HRG model, since it gathers in- 
formation, below Tc, from the whole range of possible 
imaginary chemical potentials, so that the number of 
terms actually contributing to the Fourier expansion in 
Eq. ^ can be checked with great precision: this idea 
has been followed in earlier studies limited to the 9j = 
axis [l^ . in which the presence, within errors, of a 
single Fourier contribution, corresponding to i? = 1, has 
been verified. In this respect the aim of our work is to 
extend such studies by increasing precision and by ex- 
ploring also the 9j ^ Q region. 




lm{ni) = WBjiT)S{B) cos{3B9g) 

BJ 

^ 2/3sin(2/30z) (17) 

/3>0 / 

where Wbj{T) = E^.\s.,^i'^^9,T). The average 
quark densities are always purely imaginary for imagi- 
nary chemical potentials, for that reason we shall simply 
write hq and n/ in the following, meaning implicitely that 
their imaginary part is taken. 

Predictions from the HRG model to be tested in lattice 
QCD simulations can be classified as follows: 

1) The free energy has a particularly simple form since, 
on the basis of known hadron resonances, only Wq.o, 
W^o,i: W^i,i/2i ^1,3/2 s-re different from zero in previous 
equations. That means a further strong restriction on the 
expected form of the free energy at low temperatures: a 
necessary condition for the HRG model to be valid is 
that only the few lowest terms of the Fourier expansion 
in Eq. (|4]) give contribution; 

2) Also the numerical values of the coefficients can be 
predicted from the known experimental resonance mass 
spectrum. 



FIG. 1: Normalized quark density at T = 0.9 Tc and Oj = 0. 
The solid line corresponds to the single sine fit reported in 
Table HVl 



A. r = 0.9Tc 

We start by discussing results obtained at T = 0.9 Tc- 
Let us first look at the 9i = axis: fij is zero in this 
case, while in general fig can be Fourier expanded as: 

= ^ Q sin (3Z0,) (18) 
1=1 

and the HRG model predicts contribution only from the 
lowest harmonic, 1 = 1. Indeed a simple sine term, corre- 
sponding to -B = 1, is perfectly compatible with our data, 
as showed in Fig. [1] and reported in Table IIVI A second 
term with Z = 2 is therefore not necessary, at least within 
the precision of our data, even if a two sine fit leads to 
a smaller /d-o.f. with a C2 ^ within three standard 
deviations (see again Table [IV|. As shown in Table IVl 
completely equivalent results are obtained if, instead of 
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fitting our data, we compute the coefRcients c; by explicit 
Fourier transform, 

o /•27r/3 

Q - - / sinimg)ng{0g)d9g (19) 
I" Jo 

where the integration is performed numerically by linear 
interpolation of consecutive data points. 

Next we consider data for fiq and fii obtained in the 
whole range of 9q and 0i explored, which are shown in 
Figs. [Hand [21 and try to fit them according to the expres- 
sions in Eqs. (|16p and (|17[) . considering more and more 
parameters Wbj till an acceptable value for the test 
is obtained. Fit results are reported in Table IVII a rea- 
sonable value of is obtained if a term with quantum 
numbers B = and / = 2 is allowed for, besides those 
corresponding to usual meson {B — 0, 1 — 1) and baryons 
(B = 1, / = 1/2 or 3/2). Such term does not correspond 
to any known or even possible exotic hadron [soj . but 
it is easily recognized as the first term, / = 2, neglected 
in Eq. (|12p in the Boltzmann approximation in the case 
of pions: this is actually the first correction taking into 
account quantum statistics effects for pions, i.e. the fact 
that they are bosons, and corresponds to a two-pion ex- 
change. With a pion mass as that used in our simula- 
tions, TO^ ~ 280 MeV, such term would mimic a coef- 
ficient Wo,2 ^ 0.0045, in very good agreement with the 
value obtained in our fit. Notice that terms with I > 2 
are negligible in our discretization setup, but would not 
be so, already at this temperature, in the case of physi- 
cal pion masses. As for the data at 9i = 0, allowing for 
a term with B = 2 leads to a lower value of x^, but is 
not strictly necessary, at least within the precision of our 
data. 

Our conclusion is therefore that at T = 0.9 Tc numeri- 
cal data do not contradict, within errors, the prediction 
coming from the HRG model and regarding the number 
of terms actually contributing to the free energy, apart 
from marginal evidence for a, B ^ 2 term which however 
is not strictly needed to fit data. Other deviations can 
be ascribed to the crudeness of the Boltzmann approx- 
imation for pions and are indeed well accounted for by 
the first neglected term. 

Of course if one looks at the numerical value of the 
coefficients, checking the agreement with experimental 
data is less trivial: taking into account all non-strange 
(since we are considering Nf — 2) hadron resonances 
reported in the Particle Data Book [s^, we would ex- 
pect, for instance, Wo,i = 0.457"^, which is roughly 
twice the value we have obtained (Wo,i = 0.216(2)). 
A more careful comparison is made using the unphysi- 
cal pion and p masses realized in our lattice simulations 
(m^ ~ 280 MeV and nip - 918 MeV [l^): the coeffi- 
cient becomes VFo,i ~ 0.30(2) including all resonances. 



More precisely we considered all mesons of widely accepted ex- 
istence, marked with a dot in the meson summary table. 








FIG. 2: Fit of normalized quark densities at T = 0.9 Tc, ob- 
tained from all imaginary chemical potentials explored (cross 
points), with the prediction from the HRG model (grid sur- 
face) 








FIG. 3: Fit of normalized isospin densities at T = 0.9 Tc, ob- 
tained from all imaginary chemical potentials explored (cross 
points), with the prediction from the HRG model (grid sur- 
face) 



Wo,i ~ 0.26(2) taking into account just pions and p par- 
ticles, and Wo,i ^ 0.225(15) including just pions (the 
errors here take roughly into account the uncertainties 
given for the lattice estimate of the masses in Ref. [l^), 
i.e. much closer to our numerical result or even per- 
fectly compatible in the last case. We notice that, since 
already p masses are beyond the UV scale of our lat- 
tice (a^^ ^ 700 MeV), it is perfectly reasonable that the 
contribution from higher resonances is not properly take 
into account. That also clearly shows that a compari- 
son of the numerical values of the fitted coefficients with 
the HRG model prediction is unavoidably affected by the 
systematics of the lattice discretization. 
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B. r = 0.951 Tc 

Once again we first look at results obtained for hq at 
9j — 0, which are shown in Fig. [4l In this case two 
Fourier terms, corresponding to B = 1 and B = 2, are 
necessary to fit our data. The second term is small, giving 
a contribution of the order of 5% the total signal, but our 
data are precise enough to detect it; indeed a of order 
2 is obtained if a single sine fit is tried (see Table lIVp . 

In this case the presence of the B = 2 term cannot 
be simply ascribed to a violation of the Boltzmann ap- 
proximation: assuming a mass of order 1 GeV for the 
lightest baryon, the first neglected term should lead to 
a signal a factor 10^ smaller than what we get; more- 
over it should be negative, as appropriate for a two- 
fermion exchange term. The presence in the thermal 
medium of baryon-baryon bound states, like deuterons, 
is a viable hypothesis: however assuming a mass differ- 
ence AM ~ 1 GeV between those states and the lowest 
baryon states, one would expect a suppression factor of 
the order exp(— AAf/T) ~ 10~^ at this temperature, i.e. 
much smaller than what we have obtained ^. A simpler 
explanation is that at this temperature corrections to the 
HRG model, induced by non-trivial interactions close to 
the phase transition, start to be important. 




FIG. 4: Normalized quark density at T = 0.951 Tc and Oi = 
0. The solid line corresponds to the two-sine fit reported in 
Table HVl 

That is confirmed by analyzing the complete set of 
data for fiq and fij as a function of 9q and Oj: fit results 
are reported in Table IVII Also in this case a term with 
{B,I) = (0,2) is needed, but its value comes out to be 
about twice than expected from the first term neglected 
in the Boltzmann approximation for pions. In order to 
get a reasonable value for it is necessary to introduce 



^ Notice however that also for this states lattice artifacts due to 
the low UV cutoff, ~ 700 MeV, could be important. 



also terms corresponding to i? = 2 (in agreement with 
results at 9i = 0) and terms with B — 1 and isospin up 
to / — 7/2. We interpret this again as a violation of the 
HRG model. 

Regarding the numerical values obtained for the fitted 
coefficients, we obtain for instance VFo.i ~ 0.256(2), to 
be compared with VFo,i ~ 0.24 if only pions are taken 
into account, Wo,i ~ 0.28 including pions and p mesons, 
Wq.i ^ 0.35 including all meson resonances. The same 
considerations made for T = 0.9 Tc and regarding this 
comparison also apply here. 

C. T = Tc 

Finally let us briefly discuss results obtained at T = Tc. 
Since at this temperature we stay in the confined phase 
as we switch an imaginary chemical potential, however 
small, it is still sensible to test predictions from the HRG 
model. However it is sufficient to look at results obtained 
for hq at 6j = (Fig. O to realize that violations to 
the model are important: in this case inclusion of the 
first three harmonics (B = 1, 2, 3) is necessary to obtain 
a reasonable value for (see Table IIV|) . This fact is 
confirmed by fits to the complete set of data for hq and 
hi which are reported in Table IVll the value decreases 
as more and more terms in the expansion in Eq. Q are 
added. 




0.1 0.2 0.3_ , 0.4 0.5 0.6 



FIG. 5: Normalized quark density at T = Tc and 9i — 0. 
The solid line corresponds to the three-sine fit reported in 
Table IWl 

In conclusion, within the current precision of our data, 
corrections to the HRG model are clearly detectable 
starting from T - 0.951 Tc. 

Our best fits reported in Tab. IVII which are marked in 
the x^ field by a *, provide us with a valid parametriza- 
tion of the free energy (apart from a constant term). 
We shall make use of these parametrizations in the fol- 
lowing Sections to derive generalized susceptibilities at 
fj,i = fj,q = and to study the analytic continuation of 
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the average phase of the fermionic determinant. System- 
atic effects are expected at T — Tc, where the of our 
best fit is somewhat bigger than 1. 

In order to check for systematic effects related to the 
choice of the interpolating function we have also per- 
formed polynomial fits in a limited range of chemical 
potentials: our results are reported in Table IVIII Fits 
chosen for analytic continuation are again marked by a * 
in the field. 

V. RESULTS AT T > Tc 

The range of imaginary chemical potentials available 
for analytic continuation is limited, above Tc, either by 
the presence of unphysical phase transitions related to 
center group dynamics (RW transitions) or by transitions 
corresponding to the analytic continuation of the decon- 
finement surface present at real chemical potentials. A 
full account of the high temperature phase structure in 
presence of two different imaginary chemical potentials 
will be given elsewhere [36]; in the present context we 
are just interested in the location of such transitions for 
the two temperatures explored, i.e. T — 1.048 Tc and 
T — 1.25 Tc- To that aim we have performed preliminary 
simulations on a small S'^ x 4 lattice to get a rough idea 
of the phase structure at these temperatures and thus 
delimit a safe region for analytic continuation, where to 
perform simulations on the larger 16'^ x 4 lattice. 
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FIG. 6: Polyakov loop modulus and chiral condensate at T = 
1.048 Tc along 6i = Q and 9q = Q axes. The chiral condensate 
has been divided by a factor 4 to better fit in the figure. 

As for T = 1.048 Tc, in Fig. [6] we show the behaviour 
of the modulus of the Polyakov loop and of the chiral 
condensate as a function of 6q ai 6i — Q and as a func- 
tion of 9i at 9q = Q. It is clear that along both axes a 
transition is met where the system gets back into a phase 
with confinement and chiral symmetry breaking: at those 
points the system is crossing the analytic continuation of 
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FIG. 7: Polyakov loop modulus and chiral condensate at T = 
1.25 Tc as a function of \6\ and at different constant values of 
6 1 or 6q. The chiral condensate has been divided by a factor 
2 to better fit in the figure. 

the pseudo-critical surface, present also at real chemical 
potentials. On the same symmetry grounds as for the 
deduction of general properties of the free energy in Sec- 
tion ini one expects that for small chemical potentials 
such pseudo-critical surface must be of the form 

As clear from Fig. [5] the transition happens at approxi- 
mately equal points along both axes {9q ^ Oj ^ 0.2 tt), 
i.e. A ^ B. Also the observables (Polyakov loop and 
chiral condensate) seem to be, within a good approxi- 
mation, universal functions of where 9 = {9q,9i), at 
least not too far from the origin 9q = 9j — 0. Only imag- 
inary chemical potentials strictly within the deconfined 
region can be considered for analytic continuation: Fig. [5] 

suggests us to take 16*1 < \9\max, with \9\rnax ^ 0.12 TT. 

The phase structure is less trivial at T = 1.25 Tc. 
In Fig. [7] we plot the behaviour of the modulus of the 
Polyakov loop and of the chiral condensate as a func- 
tion of 16*1 in three cases: fixed 9q = 0, fixed 9i = 
and fixed 9i = 0.257r. We observe again an approximate 
universal dependence on \9\ for relatively small values 
of this variable. Along the 9q = axis a transition is 
met, at 9j ~ 0.42 tt, which clearly belongs to the pseud- 
ocritical confinement/deconfinement surface. Along the 
9i = axis instead the system always stays in the de- 
confined phase and the Roberge- Weiss transition is met 
at 9q = 7r/3 where the system enters a different sec- 
tor, as also apparent from the behaviour of the Polyakov 
loop phase shown in Fig. [SI What happens along the 
9j — 0.257r axis is less clear: presumably there one meets 
a pseudo-critical point close to the junction between the 
Roberge- Weiss transition and the pseudo-critical decon- 
finement surface. In this context we are only interested 
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FIG. 9: Fit of normalized quark densities at T = 1.25 21:;, 
FIG. 8: Average phase of the Polyakov loop given in units of obtained from all imaginary chemical potentials in the region 
27r/3 at r = 1.25 Tc. jg| < q.SOtt (cross points), with a sixth order polynomial func- 

tion (grid surface) 



in delimiting a region safe for analytic continuation: from 
Fig. [7] it is clear that points with 16*1 < |6'|„jaj. ~ 0.3 tt 
are surely contained in that region and this has been our 
conservative choice. 

In this temperature regime we have tried to fit our 
results for hq and fii as a function of 9q , 0j accord- 
ing to polynomials derived from the general expansion 
for the free energy given in Eq. ^ and truncated to a 
given order, or according to expressions derived from a 
parametrization of the free energy given in terms of ratios 
of polynomials as in Eq. ^ . 

At T = 1.048 Tc a fourth order polynomial provides 
a good fit, while coefficients are largely indetermined if 
a sixth order polynomial is used: not enough informa- 
tion can be extracted from the limited region available 
for analytic continuation. A marginally good fit is ob- 
tained with the ratio of two second order polynomials, 
but a fourth order polynomial at the numerator seems 
preferable. 

At T = 1.25 Tc a sixth order polynomial or the ratio 
between fourth and second order polynomial are instead 
the best interpolating functions. 

A complete collection of our fit results is given in 
Table IVlIl and in Table IVTIIl Best fits chosen for analytic 
continuation are marked again by a *. Data obtained 
for T = 1.25 are shown in Figs. M a.nd\M 



VI. GENERALIZED SUSCEPTIBILITIES 

Best fits to our data provide us with a parametrization 
for the dependence of the free energy on the chemical 
potentials, from which generalized susceptibilities can be 
extracted. Results can be considered reliable as long as 
different interpolations provide consistent results. 








FIG. 10; Fit of normalized isospin densities at T = 1.25 Tc, 
obtained from all imaginary chemical potentials in the region 
\6\ < O.SOtt (cross points), with a sixth order polynomial func- 
tion (grid surface) 

In Table HI] we report results obtained for X2,o, Xi.ij 
X4.0 and X6,o (definined in Eq. ([H])) from free energy best 
fits marked by a * in the tables. Results obtained for 
X2,o, Xi,i ^^^d X4,o are reported also in Figs. [TTl [T2l and 
[13] respectively, where they are compared with analogous 
results obtained using the Taylor expansion method in 
Ref. [H. 

The following general features can be observed. Dif- 
ferent extrapolations provide always consistent results for 
X2,o and Xi,i- ^ good agreement with Taylor expansion 
results can be observed as well, apart from the T = Tc 
case. 

For X4.0 we observe a discrepancy between different 
interpolations only for T = and T = 1.048 T^; the 
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agreement with Taylor expansion is less good around Tc- 
For xg,o different extrapolations disagree or are at most 
marginally compatible in the whole range of tempera- 
tures: with the current precision of our data, we cannot 
get reliable results for sixth or higher order susceptibili- 
ties. 

In general, the comparison among different interpola- 
tion methods and with Taylor expansion results is good, 
apart from the region around Tc. This in not unex- 
pected: right above Tc the region of imaginary chemical 
potentials usable for analytic continuation is small and 
restricted by the continuation of the pseudo-critical line, 
so that poor information is available. Moreover, right at 
T — Tc we could not get best fits to the free energy de- 
pendence with a x^/d.o.f less than 1.5, therefore we do 
not have a completely satisfactory parametrization of the 
free energy for this temperature and systematic effects 
related to analytic continuation may be more important. 

We have reported in Table [J the total number of Dirac 
matrix multiplications needed in our numerical simu- 
lations at each temperature. We infer, from a rough 
estimate, that the effort for measurement purposes in 
our case (which is more or less half of the total) is ap- 
proximately two orders of magnitude larger than what 
needed (again for measurement purposes) in Ref. p^ . 
The increased effort leads to corresponding smaller er- 
rors (about one order of magnitude) only for the lowest 
susceptibilities (x20 and Xn), while for higher order sus- 
ceptibilities the Taylor expansion method seems to be 
more efficient. One has to consider, however, that our 
numerical simulations were not designed to be optimized 
for the computation of susceptibilities, and that in our 
case we obtain a complete parametrization of the free en- 
ergy dependence in terms of /ii and /i2, which is usable 
for different purposes. 

In Table imi we report also results obtained for the sus- 
ceptibilities with respect to quark and isospin chemical 
potentials and defined in Eq.®. In Fig. [T3] we show in 
particular the values of xt'o ^^'^ ^'^^ all tempera- 
tures, as obtained from polynomial fits: notice that Xo'2 
is always larger than X20 below Tc, meaning that isospin 
charge fluctuations can be excited more easily (mainly in 
the form of pions) than baryon charge fluctuations below 
Tc, while in the deconfined region the two susceptibilities 
become almost equal, as appropriate for a system made 
up mostly of quark-like degrees of freedom. 



VII. PHASE OF THE FERMIONIC 
DETERMINANT 

As we have recalled in Section II, the complex phase of 
the fermion determinant, det M[U, /x] = | det M[U, At]|e"^, 
hinders numerical simulations in presence of a real baryon 
chemical potential ^b- The problem is however milder 
in case the fluctuations of the phase (f) around zero, over 
the gauge configurations which are typical of the statis- 
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1.55(16) 


-0.0385(98) 


4.33(23) 


-69 ± 16 


1.25 


POL 


1.8470(12) 


-0.0130(9) 


1.960(20) 


0.64(23) 




RAT 


1.8473(11) 


-0.0121(7) 


1.968(16) 


2.78(25) 




\W] 


1.84(12) 


-0.0138(85) 


2.181(31) 


5.5 ± 1.7 



TABLE II: Table of different susceptibilities obtained from 
various fits. We present the values obtained from "best fits" 
of each kind of free en ergy form, together with values obtained 
by the authors of Ref. [I3] using the Taylor expansion method. 
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Fit 


q,i 
X2.0 


q,i 
X0.2 


q,i 
Xlo 


q,i 

xIa 


q,i 
X2.2 
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HRG 
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4.05(7) 


1.94(3) 




POL 


0.461(6) 


0.696(9) 


4.15(17) 
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1.76(14) 


0.951 


HRG 
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8.3(5) 


3.45(11) 




POL 
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1.589(18) 
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7.5 ±1.9 


12.8±1.9 


5.5(4) 
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3.045(21) 


3.176(15) 
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6.6(9) 


1.25 
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3.668(3) 


3.720(3) 


4.59(13) 


4.75(12) 


3.67(3) 
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3.671(3) 


3.7188(24) 


4.72(11) 


4.68(9) 


3.681(17) 



TABLE III: Table of different susceptibilities calculated with 
respect to the quark and isospin chemical potentials from the 
same best fits as for Tab. |lT] 



tical ensemble, are small: in that case efficient numerical 
methods, like reweighting, can be used. A typical mea- 
sure of the severeness of the sign problem is therefore 
given by the average of the phase factor (or some power 
of it), computed for convenience over the ensemble at 
finite isospin density. In particular in our case we can 
define: 

, ^^/2^ ^ I det Mi (/i) \ ^ Z{^l,^l) 

Z{Hg ^ = 0) 

As clear from Eq. (|20p . a way to determine (e"^/^)^ 
is to take the average of the ratio of two determinants 
over the ensemble at real isospin chemical potential: that 
is feasible but computationally demanding, especially at 
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FIG. 11: Values obtained for XiojT from various fits and 
compared with results from Ref. [l9l |. 
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FIG. 12: Values obtained for Xii/T"^ from various fits and 
compared with results from Ref. [19l ]. 



large volumes. Studying the analytic continuation of 
^gi0/2^^ at imaginary values of n, 



Z{i^, —i/i) ' 



(21) 



is an alternative: it has been shown |32l . |33| that, in 
the full QCD case, the average phase factor is analytic 
around /x^ = 0, and an efficient numerical method for the 
evaluation of the ratio of partition functions appearing in 
Eq. (pij) has been proposed in Ref. [Ill] , 

In the present context we adopt a much faster and 
cheaper approach: having measured and fitted first 
derivatives with respect to both chemical potentials, we 
have a complete knowledge, apart from constant terms, 
of the dependence of the free energy on /ii and /i2, so that 
computing the ratio in Eq. (PT|) is straightforward. Let 
us consider for instance the low temperature case, where 
we have used the HRG parametrization in Eq. (fT5)) that 



FIG. 13: Values obtained for^^o from various fits and com- 
pared with results from Ref. 
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FIG. 14: Values obtained for xlfo/T^ and Xoi/T^ from Poly- 
nomial fits. 



we rewrite: 

F = ~VT^ WBj{T)d{B) cos{3B0q 



^ 5(/3)cos(2/30/) 



(22) 



then 
Z(9„ = 



- ^ -^{F{e,=e,ei=o)-F{eg=o,ei=e)) 



(23) 



Non-zero coefficients Wb,i, apart from the constant Wo,o 
which does not enter in the computation of the average 
phase factor, have been obtained by fitting our numerical 
data. The expression can then be easily continued to real 
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chemical potentials obtaining: 

(/V3 
j;^J2^BAT)m {cos{3B9,)- 
* B,I 

^ ^(/3)cos(2/30/)) I . (24) 

/3>0 / 

The same procedure applies to other functional forms 
used in our fits: the comparison of different extrapola- 
tions to real chemical potentials based on different fitting 
functions, when available, gives a measure of the sys- 
tematic effects involved in analytic continuation. Notice 
that in the case of the HRG parametrization we can dis- 
tinguish the different contributions to the average phase 
factor, hence to the sign problem, coming from different 
particle species: this feature will be useful in our analysis. 

In Figs. [TSl and [TBI we report, as a function of 2fi/mTr, 
results obtained respectively dX T — 0.9 Tc and T — 
0.951 Tc using HRG inspired and polynomial interpola- 
tions. Where visible, the two lines reported for each ex- 
trapolation delimit the 90% confidence level region and 
give an estimate of our uncertainties: a good agreement 
between HRG inspired and polynomial extrapolations 
can be appreciated. 

It is interesting to make a direct comparison of our 
results with predictions coming from chiral perturbation 
theory (xPT). The average phase factor has been com- 
puted to one loop order of xPT in Ref. According 
to the results reported in Section VI of Ref. our 
spatial lattice size is big enough (Lsm,r ~ 6.6) to justify 
taking the thermodynamical limit at fixed T of the one 
loop xPT result, which coincides with the prediction of 
a HRG model including only pions: 

(e**/2),, = e-^°° (25) 

with 

AGo = yr3(^) J ^ cosh(2^n)-l)(26) 

This prediction (assuming in our case to^ ~ 280 MeV 
and Tc ~ 170 MeV) is reported in Figs. [I5] and [H] as 
a solid line. It is apparent that the agreement of xPT 
with the analytic continuation of our data is not satis- 
factory. In particular analytic continuation provides a 
higher value for (e"^/^)^, meaning a milder sign problem. 
To better understand the origin of this discrepancy, we 
have tried to compute the average phase factor from our 
HRG model best fit, but neglecting all contributions to 
the free energy with B ^ Q, which cannot be taken into 
account by xPT, i.e. taking only contributions from Wo,i 
and Wo, 2 in Eq- ([24]) . Results are shown in Figs. [TSl and 
[TBI in this case the agreement with xPT is almost perfect 
for T = 0.9 Tc, and acceptable for T = 0.951 T^. This 
is expected since, as we have discussed in Section flV Al 



the coefficients Wo,i and W0.2 obtained by our fits are 
compatible within errors, at T = 0.9 Tc, with those pre- 
dicted if only pions are taken into account: of course that 
may be an accident and the contribution of higher meson 
resonances should be better understood. 

Anyway, an outcome of our analysis, which is in agree- 
ment with HRG model expectations, is that contributions 
to the average phase factor coming from physical states 
with B ^ Q are significant and tend in general to make 
the sign problem less severe. 

In Fig. [T7] we report the analytic continuation of the 
average phase factor obtained at all temperatures from 
a polynomial fit: of course results reported in the fig- 
ure must be intended to be valid for chemical potentials 
bounded, below Tc, by the deconfinement critical line 
present at real chemical potentials. As expected, at fixed 
chemical potential the sign problem is much milder for 
T > Tc. This can be put again in connection with the 
fact that states with B ^ 0, which are more easily created 
above Tc, tend to mitigate the sign problem. 
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FIG. 15: The average phase factor continued from different 
interpolations and compared to 1-loop xPT results for T = 
0.9 Tc- In particular we show the 90% confidence level band 
extrapolated from our best fits to the free energy dependence. 



VIII. CONCLUSIONS 

In this paper we have studied Nf = 2 QCD thermody- 
namics, exploiting analytic continuation from two imag- 
inary chemical potentials coupled to baryon and isospin 
charges. Simulations have been performed at five tem- 
peratures around the critical value Tc ~ 170 MeV, using 
a 16'^ X 4 lattice with a standard staggered action and a 
fixed pion mass — 280 MeV. 

We have computed free energy first derivatives with re- 
spect to the chemical potentials (quark number densities) 
and interpolated them by suitable functions, in order to 



13 




2ii/m 

" 71 



FIG. 16; Same as in Fig.[T5]for T = 0.951 




FIG. 17: The average phase factor continued from polynomial 
interpolations at all explored temperatures. For each temper- 
ature we show the 90% confidence level band corresponding 
to our best fits. 

perform analytic continuation. In particular we have 
tested HRG predictions below Tc, reconstructed gener- 
alized susceptibilities at zero chemical potentials and de- 
termined the analytic continuation of the average phase 
factor. 

We have checked that HRG model predictions are in 
very good agreement with our numerical results for T — 
0.9 Tc- Small but clearly detectable deviations start to be 



visible at T = 0.95 Tc, in agreement with similar results 
reported in Ref. [13] ■ They appear, in a HRG inspired 
parametrization of the free energy, as contributions from 
unphysical states with higher values of baryon or isospin 
charges, which are of the order of a few percent at T = 
0.95 Tc and above 10 % at T ~ T^. 

Regarding the computation of generalized susceptibil- 
ities, analytic continuation gives consistent results which 
are in agreement with those obtained by the Taylor ex- 
pansion method, apart from temperatures in correspon- 
dence or right above Tc, where the range of imaginary 
chemical potentials available for analytic continuation is 
small and larger systematic effects are expected. Poor in- 
formation has been obtained for susceptibilities beyond 
sixth order. 

We have obtained consistent determinations, by ana- 
lytic continuation with different interpolating functions, 
of the average phase factor. In particular below Tc, in the 
case of HRG inspired interpolations, we have been able 
to distinguish the contribution to the average phase fac- 
tor coming from the different hadron states: results from 
analytic continuation are consistent with xPT results, 
below Tc, if one takes into account only meson contribu- 
tions. Baryons give contributions to the average phase 
factor which in general tend to make the sign problem 
less severe. The sign problem is much milder for T > Tc, 
and this can be put again in connection with the fact 
that states with B 0, which are more easily created 
above Tc, tend to mitigate the sign problem. 

Our results should be refined and could be improved in 
several respects. Simulation closer to the continuum limit 
and possibly closer to the physical quark mass spectrum 
would clarify the comparison with HRG prcdicitions, as 
well as that with xPT for the average phase factor. An 
improvement in the determination of generalized suscep- 
tibilities could be obtained by combining analytic contin- 
uation with other techniques: for instance fixing lowest 
order terms in a polynomial expansion by the Taylor ex- 
pansion method or by reweighting could lead to enhanced 
predictivity for analytic continuation. We shall continue 
our investigation along those lines in the future. 
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TABLE VI: Coefficients of HRG model fits at various temperatures. 
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TABLE VIL Coefficients of polynomial fits at various temperatures. Each line contains results of fit performed on all points in 
the circumference of radius |6|max 
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TABLE VIII: Coefficients of rational fits at various temperatures. 



